############################################################################
############################################################################
##This file produces all figures related to the Bechtel and Scheve
##conjoint example in Ratkovic and Tingley (2016).  It has 4 components:
## 1) Fit the files if the output isn't in the current directory.  Gives seed.
## 2) Produces interaction plots.
## 3) Produces reordered interaction plot so output matches with cjoint output.
## 4) Runs cjoint and produces pdf of output.
##The file can be sourced and should be run without issue.
############################################################################
############################################################################

rm(list=ls())

library(sparsereg)
library(cjoint)

load("SetupPNAS-All.RData")


##If true, prints to directory for paper
forpaper<-FALSE

############################################################################
############################################################################
##Load output if it already exists, else fit the models.
############################################################################
############################################################################

set.seed(02138)


if(file.exists("Climate-PNAS-Sparsereg.RData")) {
	load("Climate-PNAS-Sparsereg.RData")	
} else{
	g<-sparsereg(y.slim, x.slim, treat.slim,id=ID.slim, scale.type="TX",  gibbs=1000, burnin=1000, thin=30)
	save(g,file="Climate-PNAS-Sparsereg.RData")
	}
	
	
if(file.exists("Climate-PNAS-Sparsereg-EM.RData")) {
	load("Climate-PNAS-Sparsereg-EM.RData")	
} else{
	g<-sparsereg(y.slim, x.slim, treat.slim,id=ID.slim, scale.type="TX",  gibbs=1000, burnin=1000, thin=30)
	save(g.EM,file="Climate-PNAS-Sparsereg-EM.RData")	
	}

############################################################################
############################################################################
##Violin plots in paper
############################################################################
############################################################################

pdf("Figure8_PNASVolcano.pdf")
violinplot(g,columns=c("Ideology: Conservative x WhoMonitors: Your government","Ideology: Conservative x WhoMonitors: United Nations","Ideology: Conservative x WhoMonitors: Greenpeace"),newlabel=c("Conservative x Own Gov","Conservative x UN","Conservative x Greenpeace"))
dev.off()

if(forpaper){
pdf("../../../paper/figs/Figure8_PNASVolcano.pdf")
violinplot(g,columns=c("Ideology: Conservative x WhoMonitors: Your government","Ideology: Conservative x WhoMonitors: United Nations","Ideology: Conservative x WhoMonitors: Greenpeace"),newlabel=c("Conservative x Own Gov","Conservative x UN","Conservative x Greenpeace"))
dev.off()
}

############################################################################
############################################################################
##Effect plots, unordered
############################################################################
############################################################################

load("Climate-PNAS-Sparsereg.RData")


##Clean up misspellings
colnames(g$beta.mode)<-colnames(g$beta.mean)<- 
 gsub("dolar","dollar",colnames(g$beta.mode))
colnames(g$beta.mode)<- colnames(g$beta.mean)<- 
 gsub("female","Female",colnames(g$beta.mode))
colnames(g$beta.mode)<- colnames(g$beta.mean)<- 
gsub("low","Low",colnames(g$beta.mode))

##Interaction plot in paper--See below for main plot
##ggplot seems buggy; try to run the commands several times
pdf("Figure7_InteractionEffects.pdf")
plot(g,plot.one=2)
dev.off()

if(forpaper){
pdf("../../../paper/figs/Figure7_InteractionEffects.pdf")
plot(g,plot.one=2)
dev.off()
}

############################################################################
############################################################################
##Main effect plots, ordered -- 
##  reordered main effect plot to match cjoint output; requested by reviewer
############################################################################
############################################################################



load("Climate-PNAS-Sparsereg.RData")

g2<-g

##Clean up misspellings
 colnames(g2$beta.mode)<- colnames(g$beta.mode)<-  colnames(g2$beta.mean)<- colnames(g$beta.mean)<- 
 gsub("dolar","dollar",colnames(g2$beta.mode))
  colnames(g2$beta.mode)<- colnames(g$beta.mode)<-  colnames(g2$beta.mean)<- colnames(g$beta.mean)<- 
 gsub("female","Female",colnames(g2$beta.mode))
   colnames(g2$beta.mode)<- colnames(g$beta.mode)<-  colnames(g2$beta.mean)<- colnames(g$beta.mean)<- 
 gsub("low","Low",colnames(g$beta.mode))
 
 
g2$beta.mode<-g$beta.mode[,apply(g$beta.mode,2,median)!=0]
g2$beta.mean<-g$beta.mean[,apply(g$beta.mode,2,median)!=0]
g2$beta.ci<-g$beta.ci[,apply(g$beta.mode,2,median)!=0]
g2$X<-g$X[,apply(g$beta.mode,2,median)!=0]

x<-g2
main1<-"Main Effects";xlabel="Effect"
names.est <- colnames(x$beta.mode)
sum1 <- summary(x, printit = FALSE)$table
inter.count <- unlist(lapply(sapply(names.est, strsplit, 
                                    " x "), length))
length.names <- sapply(names.est, nchar)
plot.mat <- as.matrix(x$beta.mode)
num.mat <- name.mat <- plot.mat
for (i in 1:ncol(plot.mat)) {
  name.mat[, i] <- colnames(plot.mat)[i]
  num.mat[, i] <- i
}
plot.vec <- as.vector(plot.mat)
num.vec <- as.vector(num.mat)
name.vec <- as.vector(name.mat)
plotobj <- data.frame(names.est, sum1[, 1:3])
names(plotobj) <- c("names", "median50", "lo", "hi")
drop.zero <- (plotobj[, 2] != 0)
d.nz <- plotobj[drop.zero, ]
inter.nz <- inter.count[drop.zero]
lo <- hi <- median50 <- NULL
d.nz<- data.frame(lapply(d.nz,FUN=function(x) x[1:16]))


order.lev<-c(rev(c(4,1,2,3,6,5,7,8,9,10,11,12,16,14,15,13)))
d.nz$names<-factor(d.nz$names,levels=d.nz$names[order.lev])

 

g1 <- ggplot(d.nz, aes(x = names, 
                       y = median50)) + geom_errorbar(aes(ymin = lo, 
                                                          ymax = hi), width = 0.1) + coord_flip() + theme(legend.position = "none") + 
  ylab(xlabel) + theme(axis.title.y = element_blank()) + 
  geom_point(shape = "|") + ggtitle(main1) + ylim(range(d.nz[, 
                                                             -1]))+theme_bw()

pdf("Figure6_MainEffects_reord.pdf")
plot(g1)
dev.off()



if(forpaper){
pdf("../../../paper/figs/Figure6_MainEffects_reord.pdf")
plot(g1)
dev.off()
}


############################################################################
############################################################################
##Run conjoint analysis
############################################################################
############################################################################

load("pnas_cjoint.Rdata")


## Run AMCE estimator from cjoint

results <- amce(choice_cj.n ~ Cost+DistributionOfCosts+CountriesParticipating+EmissionCuts+Sanctions + WhoMonitors, data=nonmiss, cluster=TRUE, respondent.id="ID", design="uniform")


## Figure in paper
l<-  theme_bw() + theme(legend.position="none")

theme_bw_new <- function(base_size = 11, base_family = ""){
  theme_bw(base_size = base_size, base_family = base_family) %+replace%
    theme(axis.text.x = element_text(size = base_size*.9, colour = "black",  hjust = .5 , vjust=1),
          axis.text.y = element_text(size = base_size, colour = "black", hjust = 0 , vjust=.5 ), 
          
          axis.ticks = element_line(colour = "grey50"),
          axis.title.y =  element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
          plot.title = element_text(face = "bold"),
          legend.position = "none")
  
}

l<-theme_bw_new()


##Print pdfs
pdf("Figure5_cjoint.pdf")
plot(results,colors="black",plot.theme=l)
dev.off()

if(forpaper){
pdf("../../../paper/figs/Figure5_cjoint.pdf")
plot(results,colors="black",plot.theme=l)
dev.off()
}














